home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1994 August / August CD.bin / Shareware / Education / numericalmethods Folder / chap_9 / a9_s.m < prev    next >
Encoding:
Text File  |  1994-06-05  |  2.0 KB  |  79 lines  |  [MATS/MATL]
  1. echo off;
  2. % NUMERICAL METHODS: MATLAB Programs, (c) John H. Mathews 1994
  3. % To accompany the text:
  4. % NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
  5. % Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
  6. % This free software is complements of the author.
  7.  
  8. % Algorithm 9.S (Runge-Kutta Method for a System).
  9. % Section    9.7,    Runge-Kutta Methods, Page 475
  10. echo on; clc; format long; hold off; clear
  11.  
  12. % This program implements the Runge-Kutta method
  13. % for solving the initial value problem
  14.  
  15. %     Z' = F(t,Z)     with    Z(a) = Za
  16.  
  17. % where the vector notation is equivalent to:
  18.  
  19. %     x' = f(t,x,y)   with    x(a) = xa
  20. %     y' = g(t,x,y)           y(a) = ya
  21.  
  22. % Formulas for f(t,x,y) and g(t,x,y) are used to form
  23. %    the vector function F(t,Y)    which is stored in  Fn.m 
  24.  
  25. % function W = fn(t,Z)
  26. % x = Z(1);  y = Z(2);
  27. % W = [(x - x*y - x^2/10), (x*y - y - y^2/20)];
  28.  
  29. delete fn.m
  30. diary fn.m; disp('function W = fn(t,Z)');...
  31.             disp('x = Z(1);  y = Z(2);');...
  32.             disp('W = [(x - x*y - x^2/10), (x*y - y - y^2/20)];');...
  33. diary off;
  34.  
  35. fn(0,[0 0]); % Remark. fn.m and rks4 are used in Algorithm 9.S
  36. pause           % Press any key to continue.
  37.  
  38. clc;
  39. %    Place the endpoints of [a,b] in  a  and  b.
  40.  
  41. %    Place the initial value Za = Z(a) in Za.
  42.  
  43. %    Place the number of subintervals in  m.
  44.  
  45. a  = 0;
  46. b  = 15;
  47. Za = [2  1];
  48. m  = 150;
  49. [T,Z] = rks4('fn',a,b,Za,m);
  50. P = [T;Z']';
  51. points = P(1:10:length(P),:);
  52.  
  53. pause    % Press any key to see the list of solution points.
  54.  
  55. clc;, clg;
  56. X = Z(:,1);
  57. Y = Z(:,2);
  58. a = min(X)-0.3;
  59. b = max(X)+0.3;
  60. c = min(Y)-0.3;
  61. d = max(Y)+0.3;
  62. axis([0 2.1 0 1.8]);...
  63. plot(X,Y,'g');...
  64. hold on;...
  65. plot([0 2.1],[0 0],'b',[0 0],[0 1.8],'b');...
  66. xlabel('X');...
  67. ylabel('Y');...
  68. title('Runge-Kutta solution to Z` = F(t,Z)');...
  69. grid;...
  70. axis;...
  71. hold off;...
  72. shg; pause    % Press any key to continue.
  73.  
  74. Mx1 = 'Runge-Kutta solution to Z` = F(t,Z).';
  75. Mx2 = '     t(k)               x(k)               y(k)';
  76. clc,echo off,diary output,...
  77. disp(''),disp(Mx1),...
  78. disp(''),disp(Mx2),disp(points),diary off,echo on
  79.